A new test for chaos 
in deterministic systems 

By Georg a. Gottwald f and Ian Melbourne f 

We describe a new test for determining whether a given deterministic dynamical 
system is chaotic or nonchaotic. In contrast to the usual method of computing the 
maximal Lyapunov exponent, our method is applied directly to the time series data 
and does not require phase space reconstruction. Moreover, the dimension of the 
dynamical system and the form of the underlying equations is irrelevant. The input 
is the time series data and the output is or 1 depending on whether the dynamics 
is non-chaotic or chaotic. The test is universally applicable to any deterministic 
dynamical system, in particular to ordinary and partial differential equations, and 
to maps. 

Our diagnostic is the real valued function p{t) — 0(x(s)) cos(6'(s))ds where 

is an observable on the underlying dynamics x(i) and 9{t) — ct + 0(x(s))ds. The 
constant c > is fixed arbitrarily. We define the mean-square-displacement M(t) 
ioT p{t) and set K = limt^oo logM(t)/ logt. Using recent developments in ergodic 
theory, we argue that typically K = signifying nonchaotic dynamics ov K — 1 
signifying chaotic dynamics. 

Keywords: Chaos, deterministic dynamical systems, Lyapunov exponents, 
mean square displacement, Euclidean extension 

1. Introduction 

The usual test of whether a deterministic dynamical system is chaotic or nonchaotic 
is the calculation of the largest Lyapunov exponent A. A positive largest Lyapunov 
exponent indicates chaos: if A > 0, then nearby trajectories separate exponentially 
and if A < 0, then nearby trajectories stay close to each other. This approach has 
been widely used for dynamical systems whose equations are known (Abarbancl et 
al. 1993; Eckmann et al. 1986; Parker & Chua 1989). If the equations are not known 
or one wishes to examine experimental data, this approach is not directly applicable. 
However Lyapunov exponents may be estimated (Wolf et al. 1985; Sana & Sawada 
1985; Eckmann et al. 1986; Abarbanel et al. 1993) by using the embedding theory 
of Takens (1981) or by approximating the linearisation of the evolution operator. 
Nevertheless, the computation of Lyapunov exponents is greatly facilitated if the 
underlying equations are known and are low-dimensional. 

In this article, we propose a new 0-1 test for chaos which does not rely on 
knowing the underlying equations, and for which the dimension of the equations 
is irrelevant. The input is the time series data and the output is either a or a 1 
depending on whether the dynamics is nonchaotic or chaotic. Since our method is 
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applied directly to the time series data, the only difference in difficulty between 
analysing a system of partial differential equations or a low-dimensional system of 
ordinary differential equations is the effort required to generate sufficient data. (As 
with all approaches, our method is impracticable if there are extremely long tran- 
sients or once the dimension of the attractor becomes too large.) With experimental 
data, there is the additional effect of noise to be taken into consideration. Wc briefly 
discuss this important issue at at the end of this paper. However, our aim in this 
paper is to present our findings in the situation of noise-free deterministic data. 



To describe the new test for chaos, we concentrate on the continuous time case 
and denote a solution of the underlying system by x(t). The discrete time case is 
handled analogously with the obvious modifications. Consider an observable 0(x) 
of the underlying dynamics. The method is essentially independent of the actual 
form of — almost any choice of will suffice. For example if x = {xx,X2, ■ ■ ■ , Xn) 
then (^(x) = a;i is a possible and simple choice. Choose c > arbitrarily and define 



(Throughout the examples in §3 and §4 we fix c = 1.7 once and for all.) We claim 
that 

(i) p{t) is bounded if the underlying dynamics is nonchaotic and 

(ii) p{t) behaves asymptotically like a Brownian motion if the underlying dynam- 
ics is chaotic. 

The definition of p in (2.1), which involves only the observable (/)(x), highlights 
the universality of the test. The origin and nature of the data which is fed into the 
system (2.1) is irrelevant for the test, and so is the dimension of the underlying 
dynamics. 

Later on, we brieffy explain the justification behind the claims (i) and (ii). For 
the moment, we suppose that the claims are correct and show how to proceed. 

To characterise the growth of the function p{t) defined in (2.1), it is natural to 
look at the mean square displacement (MSD) of p{t), defined to be 



If the behaviour of p{t) is Brownian, i.e. the underlying dynamics is chaotic, then 
M{t) grows linearly in time; if the behaviour is bounded, i.e. the underlying dynam- 
ics is nonchaotic, then M{t) is boimdcd. Wc define K = limj^oc logM(i)/logi as 
the asymptotic growth rate of the MSD. The growth rate K can be numerically de- 
termined by means of linear regression of logM(f) versus logt (Press et al. 1992). 
(To avoid negative logarithms we actually calculate lim^^oo log(M(f) -|- l)/logt 



2. Description of the 0—1 test for chaos 




(2.1) 




(2.2) 
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which obviously does not change the slope K.) This allows for a clear distinction 
of a nonchaotic and a chaotic system as K may only take values K = ov K = 1. 
Wc have lost though the possibility of quantifying the chaos by the magnitude of 
the largest Lyapunov exponent A. 

Numerically one has to make sure that initial transients have died out so that the 
trajectories arc on (or close to) the attractor at time zero, and that the integration 
time T is long enough to ensure T ^ t. 

3. An example: the forced van der Pol oscillator 

We illustrate the 0-1 test for chaos with the help of a concrete example, the forced 
van der Pol system, 

X2 = —d{xf — l)x2 — xi + acoswt (3.1) 

which has been widely studied in nonlinear dynamics (van der Pol 1927; Guck- 

cnhcimcr & Holmes 1990). For fixed a and d, the dynamics may be chaotic or 
nonchaotic depending on the parameter w. Following Parlitz & Lauterborn (1987), 
we take a = d = 5 and let u vary from 2.457 to 2.466 in increments of 0.00001. 
Choose 0(x) = xi + X2 and c = 1.7. Wc stress that the results are independent of 
these choices for all practical purposes. As described below in §5, almost all choices 
will work. (Deliberately poor choices such as c = 0, or (/> = 7 or ^ = t obviously 
fail; sensible choices that fail are virtually impossible to find.) 

In figure 1 we show a plot of K versus w. The periodic windows are clearly seen. 
As a comparison we show in figure 2 the largest Lyapunov exponent A versus w. 
Since the onset of chaos docs not occur until after ui = 2.462 we display the results 
only for the range 2.462 < a; < 2.466 in figures 1 and 2. (Both methods easily 
indicate regular dynamics for 2.457 < u) < 2.462.) 

Naturally we do not obtain the values K = and K = 1 exactly - the mathe- 
matical results that underpin our method predict these values in the limit of infi- 
nite integration time. (The same caveat applies equally to the Lyapunov exponent 
method.) In producing the data for figures 1 and 2, we allowed for a transient of 
200,000 units of time and then integrated up to time T = 2, 000, 000. As can be seen 
in figure 3, for most of the 400 data points in the range of w, we obtain > 0.8 or 
K < 0.01. 

Next, we carry out the 0-1 test for the forced van der Pol system in the situation 
of a more limited quantity of data. The results are shown in figure 4 for 2.463 < 
u! < 2.465. Wc again allow for a transient time 200,000 but then integrate only for 
T = 50, 000. The transitions between chaotic dynamics and periodic windows are 
almost as clear with T = 50, 000 as they are with T = 2, 000, 000 even though the 
convergence of to or 1 is better with T = 2, 000, 000. 
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Figure 1. Asymptotic growth rate K of the mean square displacement (2.2) versus to for 
the van der Pol system (3.1) determined by a numerical siirmlation of the skew product 
system (3.1) and (2.1) with a = d = 5, c = 1.7, <^(x) = x\ X2 and uj varying from 2.462 
to 2.466. The integration interval is T = 2, 000, 000 after an initial transient of 200, 000 
units of time. 
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Figure 2. Largest Lyapunov exponent A versus uj for the van der Pol system (3.1) with 
o = d = 5 and uj varying from 2.462 to 2.466 (of Parlitz & Lauterborn 1987). The 
integration interval is T = 2, 000, 000 after an initial transient of 200, 000 units of time. 



Article submitted to Royal Society 



A new test for chaos 



5 



* + + 
*+ *+ +++ 



-■•»#»>»+""^""+^*"1 



^^^^ ^ 



Figure 3. Asymptotic growth rate K versus u> for the van der Pol system (3.1) as in 
figure 1 with T = 2, 000, 000. The horizontal Unes represent K = 0.01 and K = 0.8. 
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Figure 4. Asymptotic growth rate K versus uj varying from 2.463 to 2.465 for the van der 
Pol system as in figures 1 and 3 but with integration interval T = 50, 000 after an initial 
transient of 200, 000 units of time. The horizontal lines represent K = 0.01 and K = 0.8 . 
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4. Further examples 

To test the method on high-dimensional systems we investigated the driven and 
damped Kortweg-de Vries (KdV) equation (Kawahara & Toh 1988) 

Ut + uu^ + /3u^xx + au^x + vu^xxx = 0, (4.1) 

on the interval [0, 40] with periodic boundary conditions. This partial differential 
equation has non-chaotic solutions if the dispersion /3 is large and exhibits spatio- 
temporal chaos for sufficiently small /3. Note that equation (4.1) reduces to the 
KdV equation when a = u = 0, and reduces to the Kuramoto-Sivashinksy equation 
when /3 = 0. 

We fix a = 2, ly = 0.1 and vary /3. For /3 = 0, it is expected that the dynamics 
of the Kuramoto-Sivashinksy equation are chaotic for these parameter values. As 
an observable we used (/)(u(a;, t)) — u{xo, t) where xq is an arbitrarily fixed position, 
and we iterated until time T = 35, 000. The 0-1 test confirms that the dynamics 
is chaotic at /3 = (with K = 0.939). Also, wc obtain K = 0.989 at /3 = 0.1 and 
K = 0.034 at /? = 4, indicating chaotic and regular dynamics respectively at these 
two parameter values. 

Finally, for discrete dynamical systems, we tried out the test with an ecological 
model whose chaotic component is coupled with strong periodicity (Brahona & 
Poon 1996; Gazelles & Ferriere 1992). The model 

Xk+i = 118j/feexp(-0.001(a;fc -l-t/fe)) 

Uk+i = 0.2 Xk exp(-0.07(a;fe +yk))+ 0.8 yk exp(-0.05(0.5 Xk + Vk)) 

has a non-connected chaotic attractor consisting of seven connected components. 
Our tost yields K = 1.023 with only 10, 000 data points and clearly shows that the 
dynamics is chaotic. 

5. Justification of the 0—1 test for chaos 

The function p{t) can be viewed as a component of the solution to the skew product 
system 

e = c + </>(x(t)) 

p = (j){x{t)) cose (5.1) 
q = ^(x(f)) sin^ 

driven by the dynamics of the observable (/)(x(t)). Here {6. p, q) represent coordinates 
on the Euclidean group E(2) of rotations 9 and translations {p, q) in the plane. 

We note that inspection of the dynamics of the (p, <7)-trajectories of the group 
extension provides very quickly (for small T) a simple visual test of whether the 
underlying dynamics is chaotic or nonchaotic as can be seen from figure 5. 
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Figure 5. The dynamics of the translation components (p, g) of the E(2)-extension (5.1) 
for the van der Pol system (3.1) with a = d = 5, c = 1.7, </>(x) — xi + X2- These plots 
were obtained by integrating for T ~ 1400 (with timestep 0.01). In (a), an unbounded 

trajectory is shown corresponding to chaotic dynamics at a; = 2.46550. In (b), a bounded 
trajectory is shown corresponding to regular dynamics at w = 2.46551. 
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In Nicol et al. (2001) it has been shown that typicaUy the dynamics on the group 
extension is sublinear and furthermore that typically the dynamics is bounded if the 
underlying dynamics is nonchaotic and unbounded (but sublinear) if the underlying 
dynamics is chaotic. Moreover, the p and q components each behave like a Brownian 
motion on the line if the chaotic attractor is uniformly hyperbolic (Field et al. 
2002). A nondcgcncracy result of Nicol et al. (2001) ensures that the variance of 
the Brownian motion is nonzero for almost all choices of observable (j). Recent work 
of Melbourne & Nicol (2002) indicates that these statements remain valid for large 
classes of nonuniformly hyperbolic systems, such as Hcnon-likc attractors. 

There is a weaker condition on the 'chaoticity' of X that guarantees the desired 
growth rate = 1 for the mean square displacement (2.2): namely that the autocor- 
relation function of (/)(x) cos(0) decays at a rate that is better than quadratic. More 
precisely, let x(t) and 9{t) denote solutions to the skew product equations (5.1) with 
initial conditions xq and respectively. If there are constants C > and a > 2 
such that 



for all t > 0, then K = I as desired (Biktashev & Holden 1998; Ashwin et al. 2001; 
Field et al. 2002). (Again, results of Nicol et al. (2001); Field et al. (2002); and 
Melbourne & Nicol (2002) ensure that the appropriate nondegeneracy condition 
holds for almost all choices of </>.) There is a vast literature on proving decay of 
correlations (Baladi 1999) and this has been generalised to the equivariant setting 
for discrete time by Field et al. (2002) and Melbourne & Nicol (2002) and for 
continuous time by Melbourne & Torok (2002). It follows from these references 
that K = 1, ioY large classes of chaotic dynamical systems. 

One might ask why it is not better to work, instead of the E(2)-extension, with 
the simpler M-extension 



In principle, p{t) can again be used as a detector for chaos. However, by the ergodic 
theorem p{t) will typically grow linearly with rate equal to the space average of 
(j). This would lead to K = 2 irrespective of whether the dynamics is regular or 
chaotic. Hence, it is necessary to subtract off the linear term of p{t) in order to 
observe the bounded/diffusive growth that distinguishes between regular/chaotic 
dynamics. Subtracting this linear term is a highly nontrivial numerical obstruc- 
tion. The inclusion of the 9 variable in the definition (2.1) of p{t) and in the skew 
product (5.1) kills off the linear term. 



We have established a simple, inexpensive, and novel 1 test for chaos. The com- 
putational effort is of low cost, both in terms of programming efforts and in terms 
of actual computation time. The test is a binary test in the sense that it can only 
distinguish between nonchaotic and chaotic dynamics. This distinction is extremely 
clear by means of the diagnostic variable K which has values either close to or 
close to 1. The most powerful aspect of our method is that it is independent of the 
nature of the vector field (or data) under consideration. In particular the equations 




(j){Ti.{t)) cos{e{t))(l){xo) coseodxod6'o < Ct' 



.—a 



p = .^(x(i). 



6. Discussion 
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of the underlying dynamical system do not need to be known, and there is no prac- 
tical restriction on the dimension of the underlying vector field. In addition, our 
method applies to the observable (f>{x{t)) rather than the full trajectory x(t). 

Related ideas (though not with the aim to detect chaos) have been used for 
PDE's in the context of demonstrating hypermeander of spirals in excitable me- 
dia (Biktashcv & Holdcn 1998) where the spiral tip appears to undergo a planar 
Brownian motion (see also Ashwin et al. 2001). We note also the work of CouUet 
& Emilsson (1996) who studied the dynamics of Ising walls on a line, where the 
motion is the superposition of a linear drift and Brownian motion. (This is an ex- 
ample of an R-extension which we mentioned briefly in §5. As we pointed out then, 
the linear drift is an obstruction to using an M-extension to detect chaos.) 

From a purely computational point of view, the method presented here has a 
number of advantages over the conventional methods using Lyapunov exponents. 
At a more technical level, we note that the computation of Lyapunov exponents 
can be thought of abstractly as the study of the GL(n)-extension 

A = (df)x(t)A 

where GL(n) is the space of n x n matrices, A G GL(n), and n is the size of the 
system. Note that the extension involves additional equations and is defined us- 
ing the linearisation of the dynamical system. To compute the dominant exponent, 
it is still necessary to add n additional equations, again governed by the linearised 
system. In contrast, our method requires the addition of two equations. 

In this paper, we have concentrated primarily on the idealised situation where 
there is an in principle unlimited amount of noise-free data. However, in §3 we also 
demonstrated the effectiveness of our method for limited data sets. An issue which 
will be pursued in further work is the effect of noise which is inevitably present in 
all experimental time series. Preliminary results show that our test can cope with 
small amounts of noise without difficulty. A careful study of this capability, and 
comparison with other methods, is presently in progress. 

We are grateful to Philip Aston, Charlie Macaskill and Trevor Sweeting for helpful dis- 
cussions and suggestions. The research of GG was supported in part by the European 
Commission funding for the Research Training Network "Mechanics and Symmetry in 
Europe" (MASIE). 
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